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DOUBLE TIME LAG COMBUSTION INSTABILITY MODEL 


FOR BIPROPELLANT ROCKET ENGINES 


INTRODUCTION 


Since it was first observed in the early 19^*0' s, low frequency combustion 
instability or chugging in liquid rocket engines has been the subject of many 
analyses. Von Karman (l) was the first to propose that the phenomenon is due 
to a combustion time delay between the instant of propellant injection and 
subsequent conversion into combustion products. In 1950 Gunder and Friant (2) 
presented an analysis in which this combustion delay was the essential feature 
but also included the inertia of the liquid in the feed system. This was 
followed about one year later with an analysis by Summerfield (l) which also 
incorporated the combustion delay and feed system inertia. Gunder and Friant 
treated both monopropellant and bipropcllant rocket systems with a common 
combustion delay for the latter case. Summerfield treated only the monopro- 
pellant case. In both of these analyses the authors showed that instability 
is not possible if the pressure drop across the injector is greater than 
one-half the chamber pressure. Crocco and Cheng (3) later refined the model 
for the monopropellant case by assuming a time varying combustion delay which 
for simplicity was correlated to chamber pressure. Their analysis considered 
feed system inertance and capacitance as well as resistance. Hurrell (U) 
introduced the concept of an injection velocity-dependent combustion delay and 
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its effect on the neutral stability boundaries. More recently Wenzel and 
Szuch (5) conducted an analysis of the bipropellant case by allowing different 
vaporization rates for the two propellt its . However, feed system inert anc*' 
and capacitance were not considered. An interesting conclusion from this work 
was that in some cases decreasing the ratio of the injection pressure drop to 
chamber pressure results in a transition from unstable to stable operation. This 
cannot be predicted from single combustion delay analyses. 

All of these analyses have served to establish the current knowledge and 
understanding of low frequency combustion instability and have guided the pre- 
vention and/or elimination of the phenomenon in past rocket engine development 
work. However, in analyzing current high chamber pressu e, bipropellant rocket 
engines with propellants possessing distinctly different vaporization rates 
(time lags), these analyses have certain shortcomings. First, the analysis of 
a bipropellant rocket engine with a monopropellant model applied individually 
to each propellant system yields questionable results because it neglects the 
influence of the other system on the overall stability. Second, the use of a 
bipropellant model with a common combustion delay to represent propellants with 
distinctly different vaporization rates is unrealistic. And third, in view of 
the fact that feed system inertance and capacitance play an important role along 
with injection pressure drop in determining the stability of the overall system, 
these factors should also be included in the model. It should be possible, at 
lease in some cases, to stabilize the combustion with less impact to the overall 
rocket system by optimizing these parameters rather than adjusting only the 
pressure drop. Reliance on an injection pressure drop greater than one-half the 



chamber pressure to guarantee st: Mlity is unrealistic from both the standpoint 
of high pump discharge pressures that would be required for today's engines and 
from the results of the analysis or Wenzel and Szuch. 

This report advances a bipropellant stability model in vhich feed system 
inertance and capacitance are treated along with injection pressure drop and 
distinctly different propellant time lags. The model is essentially an 
extension of Crocco's and Cheng's moconropellant model to the bipropellant 
case assuming that the feed system xnertance and capacitance along with the 
resistance are located at the injector. The neutral stability boundaries are 
computed in terms of these parameters to demonstrate the interaction among 


them. 



Combustion Chamber and Associated Equations 


Crocco'a Ingenious derivation of the combustion chamber equation is 
based only on a limited number of fundamental assumptions and definitions. 
First, like Summerfield, he postulates that the chamber pressure p is a 
function of time t, even when the combustion process is a steady one. Next, 
he suggests that the rate of combustion processes in the chamber is a function 
of several variables, the two more prominent of which are pressure and temp- 
erature; the rest of the relevant variables are lumped into a single group, 

Z. Thus, without identifying the nature of the combustion process, he states 
neatly. 

Process Rate ■ ftp,T,Z) ■ f(p,T,2)[l ♦ p' (1) 


where the bars on the variables indicate the steady-state values of these 
variables, and p' ■ p - p is the small perturbation of the chamber pressure, 
and H.D.T. denotes the higher derivative terms, such as 



3f dT ^ 3f dZ. 
7T 3p 3T Hp 3 


P 


p. T = i, z 


z. 


In so doing, Crocco singles out the predominant effects of pressure on 
the process rate f, and disregards those of the others. As it will be 
shown later that, out of this perfectly general and vague definition of f, 
Crocco was able to lay the foundation for the formulation of the relationship 
between the burning and the injection rates of the propellant. 

From the definition of f just introduced, it is easily obtained by 


transposition 
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n , P'tH.D.T.) , f(p,T,Z) - f(p,T,Z) C2) 

f(p.T.Z) f(p,T,Z) 

where n is simply the percent change of process rate, with respect to the 
steady state process rate, due to a small perturbation in pressure, p'. 

Crocco and Cheng call this % change the interaction index. This index seems 
to reflect, to some degree, the design of combustion chamber. 

The next innovation that Crocco introduced at the outset of his in- 
vestigation in combustion instability is the concept of a time lag, x , 
which is the total time elapsed between the instant when n = 0 and the instant 
when n = 1. To simplify the ensuing analyses, he assumes further that during 
a certain portion of this time lag, i.e. x^, the interaction index n is zero 
and during the rest, i.e. x*x-T^,n*l. Thus, x^ is the insensitive 
part and x is the sensitive part of the total time lag x . 

With these two quantities defined, the following statement concerning 
the energy E contained in a certain element of propellant as it transforms 
from liquid state into gaseous products of combustinn can be made 



f(t') dt' * E 

a 


( 3 ) 


Note that the lower limit of the integral is the instant when the process 
rates begin to be affected by the combustion processes. Note, too, that 
this same level of energy E would have been reached also, if the processes 
of combustion had been steady. Thus, it is equally valid 
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) dt 1 =« E 

a 


i 


* 

? 


(4) 


where the schematic variation of the rate function f with respect to the 
time lags t and 7 are shown in Figure ft. 



Time 


t- T t t 

J f C t ') dt ' + / f(t’) dt ' ■ E a ■ r 7(t») dt' (5) 

t-T t-7 t-7 

The first integral on the left-hand side of (5) can be rewritten approximately, 


t-T 


/ f(t') dt* »[(t - 7 ) - (t - t)] 7(t-7) ■ (t - f) T( t - 7 ) (6) 

* t-T 


Transposing, 
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t-T 

t-x = — f f(t ) dt 

f(t-i) J 

t-T 


But, from (5), 
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t-T t 

/ f(t’) dt ' = -/ [f(t') - f(t')] dt' 

t-T t-T 

hence, (7) becomes 

tLU.:.ll£l dt ', . f - if.*) dt ' 

t-; ) 

Using the definition of n from (2), we have 

t 

n C 1 1 * 

t _ t = _ _ / p (t ) dt (8) 

P t-i 

Differentiation with respect to t of t-t under the integral sign yields then 

■ -= tp'(t) - p'(t-f)] (9) 

P 

This equation only portrays a part of the combustion chamber phenomenon. 

The second part of the derivation of the chamber equation is con- 
cerned with the mass balance in the combustion chamber. It begins with 
the premise that the mass of propellant injected equals the mass burned 
at an instant t. later. 

6m^(t) dt ■ 6m^(t+T t ) (dt + dt) (10) 

where the subscripts i and b denote injection and burning, respectively. 

The time lag T t in (10) is approximately equal to and will be later 
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replaced by its counterpart x t in the steady state process. Now let 

t+T t ■ T or t»T - x^. 

Eq. (10) becomes 

6m i (T-T t ) (dT-dx) = 6m b (T) dT 

or 

a- |f) 

Renaming T as t. as the starting instant of the process is immaterial, we 
have 

6m b (t) * &h.(t-i t ) (1 - ~) (11) 

Eq. (1) is now modified slightly, with the understanding that in a steady- 
state process fiiL ■ <$m b , 

6m fa (t) - 6in b - fim^t-x^ - 6m i - ~ Sm^t-x) (12) 

If the injection rate is constant, i.e. 6riu (t) * fiiL ■ 6m b , it follows 
then 

“bW - «*b (1 - 3F>' 

The next phase of investigation deals with the dynamics of the com- 
bustion chamber. For a non-steady process we can write 

* p (t) = * e (t) ♦ Mjt) (13) 

where m^t) is the rate of generation of combustion products, m 0 (t) is 

the rate of ejection of gases through the nozzle, and M g (t) is the mass 

of gases accumulated inside the chamber. Since M (t) is proportional to 

© 

the chamber pressure, the rate of accumulation of M^(t) is 



9- 


-3t M g")- S g 3F 

P 

where M is the steady-state value of M (t) . 

S 6 

Considering the total amount of the combustion product from t 

the current instant t c t, we have 

t t-T. 

/ i » r ^ i i 

m b (t ) dt = J m i (t ) dt 

o o 

Differentiation of (14) with respect to t gives 

A b (t) = (1 - |£) m.(t-T t ) 

Substitution of (9) into (15) yields 

m b (t) * { 1 + n [2J£X - l & L *) ]jm i (t-T t ) 

or, with a slight modification, 

itL i . ! , , .i p'w-w . lmi Ct - T ' ) ' w 

a v 


P J m 


Denoting 


m h (t)-in m.(t-T )-m 

V‘> . ■ 1 

m 


T 

m 


IK (t-r t ) and P . (*1~P .. 
1 P 


we have, 

y b (t) » n [+(t) - <Kt-i)] + ^(t-^) + n[<Kt) -♦(t-t)]y i ( 
or, approximately. 


P b ctj - ^(t-i t ) + n U(t) - ♦(t-T )] 


■ 0 to 

(14) 

( 15 ) 

— & 

- 4(t) 

t-T t ) 

(16) 
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Introducing the gas residence time 

V M ^i > (17) 

as the time an average element of products of combustion will remain in 
side the chamber in a steady- state operation before it emerges from the 
nozzle, Eq. (13) becomes 

3 I + V (Z) = P b (z) (18) 

where z - t /0 . (19) 

s 

The ejection rate can be calculated from the steady-state nozzle transfer 
function [ 6 ] as 


y e (z)as<Kz) (20) 

Substitution of (16) into (18) gives the equation of combustion chamber 
dynamics. 

+ 4 »(z) = ^(z-^) + n [*(*) - <Kz-t)] (21) 

where t has been approximated by the steady-state value t and has also 
been non-dimensional ized by the use of 0^. 

Tht term ^(z-t^) in (21) can not be determined without an examination 
of the mechanics of the feed system of the rocket engine. 

Obviously, if the injection rate is independent of the pressure 
oscillations in the chamber, ( 21 ) reduces to its simplest form. 


+ (1-n) 4> ( z) = - n 4 >(z-t) 


( 21 ) 
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Derivation of Feed System Equation 

A schematic feed system for a monopropellant rocket engine as shown as 
in Fig. 2 p can be mathematically represented by several simple component 
equations, each of which portrays a specific portion of the operation. 



Pump: 

*0 *0 

HI 1H 

_ - _ n 0 0 

(22) 


?o 

*~ L/ 

in 

0 

Mne 0-1: 

V A 1 

= p o X dt 

(23) 


SX 

1 

0 

CX 

. din 

= Zi 0 

A dt 

(24) 

Line 1-2: 

h-* 

1 

K) 


(25) 

Injector plate: 

in 2 

(26) 


(No feedback control) 


(27) 
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w’nere D and X are the proportional constat o' the pump and the equivalent 
spring constant of the feed line, respectively. X is closely related to the 

I 

feed line capacitance C^. 

H.S. Tsien [3] obtained a differential equation from (22) - (27) which 
relates the feed system dynamics with the combustion chamber phenomenon. 

f [♦ ♦ (P ♦ $ ♦ JEy 0] 

+ [l + D (P + i)] y. ♦ [DE(P ♦ i) ♦ j] 0 

i d 2 y. d V; 

+ [DJE(1 - y)(P ♦ j) ♦ JF.y] jpA- * J 2 Ey(l-y)^pr=- - 0 (28) 

where 

P * ■=• — £• (pressure drop parameter) 

Ap 

2App X 

E = — - (Elasticity parameter) 

m e g 


SL& 


2ApAt' 


(Inertia parameter) 


g 


Simplification of (28) is possible for specific vases: 
I. Constant feed pressure, D = 0 


dy. 


d 2 v 


pc* ♦ JEy p-} + ^ + J ar- + JE y 


d 2 p^ 

dP~ 


(29), 


H. Line elasticity or line capacitance concentrated at injector plate, y=l 


2 dp. d 2 p. 

P(* + p) + ^ ♦ J + ^ -0 


(29) 


II 
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III. Constant feed rate, D -*■ « 


PE 


d<j> 

dz 


du. 


* \ * E dT * (1 -”’ JB 


d 2 y . 

_i 

dz 7 " - 


(29) 


III 


IV. Line elasticity or Line capacitance concentrated at the tank end 
of line, y = 0 . 

P[* + DE(P + I) M ] + [1 + D(P + i)]y. + [DE(P + j) * J] ^ 


1 d2y i 

+ OJE (P + |) ^ = ° 


(29) 


IV 


For each special case (four sample cases are shown above) the feed 
system equation F(y^,$ ) * 0 (29) must be solved simultaneously with the 
chamber equation C($,ik) * 0 ( 21 ). if one is interested in the stability 
aspects of the problem, only the characteristic equation is of importance. 
Some sample analyses leading to this chamber equation will be given below. 

Case I. Define a differential operator 

® « 1 + J 37 ♦ JEy ♦ J 2 Ey (1 - y) - 3 - . (30) 

Then (29)^ becomes for the instant of z-t^ 

p [I * ^ 3?%-;) * ® ■ 0 < 31 > 

Applying the differential operator J) defined by (30) into (22) we have 

£> Bl - n) 4 >(z) - n<p (z-i) + * 3 y.(z-i t ) (32) 

Substitution of (32) into (31) yields a differential equation in 4 ) only. 



1 
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[1 + j jp* JEy ♦ J 2 Ey(i - y) [(!-«)♦(*) 


+ Hz ~ n ^ Z ” 


•f 


p[l ♦ JEy <Kz -*) 


0 


or 


{(1-n) ♦ [i ♦ (l-n)j] ^ + J [1 + (l-n)Ey] ^ JEy [(1-Wl-y) + l] 


+ J 2 Ey (i-y) ] 


<Kz) 


+ [P - n - nJ ♦ (P-r)JEy ^ - nj 2 dy(l-y) ^y] *(z-*) - 0 (33) 


Equation (33) can be written symbolically 


4 [♦(*) ] + 4 [«Kz-t) ] e 0 (34) 

where 

Lj = (1-n) + [1 ♦ (l-n)J] ^ * J [1 + (1-n) Ey] ^ 

♦ JEy [(1-n) J (1-y) + l] ♦ J 2 Ey (1-y) , (35) 


and 


d 2 


L 2 a (P-n) - nJ + (P-n) JEy - n J 2 Ey (1-y) 


Now we relate the function with a retarded variable (z-t) with +(z) 
by the use of Laplace transform, namely, 

CD 

*(s) - J e" sz *(z) dz 

o 


and 
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J e' sz (Kz-t) dz * e‘ ST »(s) (36) 

o 

Thus, 

j(l-n) ♦ [l + (l-n)J]s + J [l + (l-n)Ey]s 2 ♦ JEy [l + (l-n)J(l-y)]s 3 

+ J 2 Ey*,l-y)s 4 j #(s) - [l + (l-n)J] $(o) 

- [1 ♦ (l-n)Ey][s<Ko) * y 1 «] 

-JEy[l ♦ (l-n)J(l-y)][s 2 iJ>(o) ♦ s ^ | z „ 0 + | Z .J 

-J 2 Ey(l-y)[s 3 <fr(o) ♦ s 2 | 2=0 ♦ s £z*[ m0 + 3T3*| Z . 0 ] 

♦ e' ST J[(P-n) - nJs + (P-n)JEy s 2 - nJ 2 Ey(l-y)s 3 ] *(s) 

♦ nJ*(o) - (P-n)JEy [s*(o) + gf | Z . Q ♦ 

♦ nJ 2 Ey(l-y) [s 2 ^(o) + s g£| lmQ * g-£| z « 0 ]]“ 0 < 37) 

Solving for $(s), we have 

1(3) - 0 1 (s)/0 2 (s) (38) 

where 

^(s) - [l J (l-n)jJ <^(0) + [l*(l-n)Ey][sf(0)+f (0)] +JE^[l+(l-n)J(l-y)] . 

.[^(0)*#f'(0H"(0)] ♦ J*Ry(l-y) [s^(0)48 2 ^'(0)*8f'(0)^'“(0)] 

- e' 8t { n jf(0) - (P-n)JEy [af(0W (<>)]♦ nJ*Ey(l-y)[s 2 *(0)4sf • (0W'(0)]} 
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0 2 (s) ■ (l-n)+ [l+(l-n)j] s ♦ j[l+(l-n)Ey] s 2 
♦ JEy^l+(l-n)J(l--y)j 3 ^ + J^E^(l-y)s^ 

- e“ 3 * [(P-n)-nJs+(P-n}JEys 2 - nJ^Ey(l-y)s^ ] 

Eq (38) can now be inverted. The inversion process involves the 
evaluation of the residues at the poles inside a contour which encloses 
all the poles of the integrand. Since there are no poles other than those 
introduced by the vanishing of the denominator, we only need to set it to 
zero and seek its roots, either real or complex. Hence 

(l-n)+[l+(l-n) j] s+J [l+(l-n)Ey] s 2 +JEy [l *(l-n) J(l-y) ] s 3 +J 2 Ey(l-y)s 4 

= e ST [(P-n) - nJs + (P-n) JEys 2 - nJ 2 Ey (l-y)s 3 ] (39) 

If y *■ L, such that the .line capacitance is concentrated at the injector plate 
end of the line, (39) further reduces to 

(1-n) + [l + (l-n)j] s t j[l t (l-n)Ejs 2 + JES 3 -is 

P - n - nJS ♦ (P-n) JEs 2 ° e ( 4 °) 

For given values of J,E,P,N and 7 , the roots (either real or complex) 
can be solved from (40) . Since these roots will give rise in the solution 
4>(z) to such terms as 

i * 0,1,2, . . .« 


c. e 
i 


(41) 
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We are certain that a real positive root or a positive real part of a 
complex root will cause <p (z) to increase without bound as z increases. In 
other words, the existence of such a root signifies instability, whereas the 
existence of a negative real root or a negative real part of a complex 
root indicates stability. 

Case HI. Constant feed rate [>*“ 

Define ft = 1 ♦ E jL + (1-y) JE ^ (42) 

then (29)111 becomes for the instant of z - 7^ 

PE j- 2 - 4(z - t) + J»y.(z-7 t ) - 0 (43) 

But$iJ-(z - 7 ) is, from (21), 

X t 

JH'iC* - ~ t ) ■ ?C(1 - r. ♦ -Jj-MOO ♦ n *(?. - t)] (44) 

Hence 

PE jjj. +(z - 7 ) ♦ [1 ♦ E ^ + (1 - y) JE ~r][(l - n ♦ ^) *(z) + n4(z - 7 )] - 0 
Simplifying, 

(1 - n + ^)[1 + E - y) JE 

♦ {PE + n(i ♦ B jj* ♦ (1 - y) JB - t) ■ 0 (45) 

The counterpart of Eq. (39) for this case is 

(1 + Bs ♦ (1 - y) JEs 2 )[l + s - n + n e" TS ] + PEs e" TS - 0 (46) 

Eq. (46) can be solved for s for s real or for s, complex. For real s, we 
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transpose (46) as 


^-ts (1 - n + s) [1 + Es + (1 - y)JEs 2 } 
e = PEs + n[l + Es + (1 - y)JEs2) 


(47) 


Graphical procedure must be used to construct curves for the left and right 
sides of (47) for various values of the parameters n,E,y,J, and P. The 
intersections of these curves locate the roots. Again, positive roots indi- 
cate stability, otherwise, instability. 

For s complex, we substitute in (46) 

s = o + iw 

resulting in 
e" OT (cos 7w - i sin 7u>) 

[(1 - n + o)Aj - u^] + i[uiA^ + (1 - n + 0 ^ 2 ] 

" (PEa + nA x ) + i(PEw + nA ? ) 

whei-e Aj • 1 + Eo + (1 - y)JE(a 2 - u) 2 ), 

A^ 0 Eut + 2(1 - y)JEaa>. (49) 


Equate the real and imaginary parts in (48) , we obtain 


-or _ (PEa ♦ nA^) [(l-ns-ajAj-tjAgj-*- (Pfiui + (l-n+a)A 2 ] 

(PEa + •» (PEai + nA 2 ) 2 


( 50 ) 


e aX sin tea 


(PEa + nA]^) [(l-n+a)A 2 + oAjj - (PEuj + nA 2 )ftuA 2 + (l-n+a)Aj 

2 2 
(PEa + nA-^) + (PEw + nA 2 ) 


Further simplification yields 
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(PEa + nA ) cuA +(l-n+a)A 2 - (PEw + nA ) (l-n+a)A,- gjA 2 

- tan to) * ■■ - - ■' 1 — 1 

-(PEa + nA 1 ) <‘> A 2 ~( 1 “ n+a ) A i + (PEw * nA 2 ) (l-n+a)A 2 -<i>A, ^1) 


z - [(1 - n + a) 2 + u) 2 ][A 2 + A 2 ] 
(PEa + nAj) 2 + (PEid + nA 2 ) 2 


(52) 


The case of neutral stability is characterized by a * 0. Setting 
a = 0 in (51) and (52) we obtain 


nA?[u>A° + (1 - n)A°] - (PEw + nA")[(l - n)A? - ojAH 
-tan Tto * i , 

-nAj[u>A° - (1 - n)Aj] + (PE» + nA°)[(l - n)A° - cuAj] 


(53) 


[(1 - n) 2 + a) 2 ] C (A ) 2 + (A.) 2 ] 

1 » i 1 

(nAj) 2 + (PEw + nA|) 2 


(54) 


where Aj = 1 - J'Ew 2 
A 2 ■ Euj 


J' = (l - y)J. 


Simplifying, (53) , (54) becomes 


(1 - n) ♦ (P + n) [;—+ 

-tan Tu) * — ■%— - 

<*> ♦ (1 - (P ♦ n)[£ ♦ 

and 

(P + n) 2 - n 2 ■ [w 2 + (1 - n) 2 - n 2 ](^4 2 +1) 

y* jw ■ fer . 


(55) 


where 


(56) 
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THE PRESENT MODEL 

The low frequency combustion stability model described here is an 
extension of Crocco and Cheng's (4) analysis of the monopropellant system 
to the bipropellant case. The model considers the injector capacitance and 
inertance as well as the resistance of both propellant systems and allows 
for separate and distinct time lags for each propellant. Previous bi- 
propellant models do not consider the effects of injector capacitance and 
inertance. The time lags are the time intervals between the fuel and oxidizer 
injection and the assumed sudden conversion to exhaust products; they include 
all the physical and chemical processes in the conversion such as heating, 
vaporization, mixing, and reaction. 

ANALYSIS 

The monopropellant, single time lag model of Crocco and Cheng is modified 
to accommodate the bipropellant case by adding a term accounting for the 
second propellant to the equation governing the dynamics of the gas flow in 
the combustion chamber and adding a new equation representing the dynamics of 
the second feed system. The modified equation for the chamber dynamics in 
dimensionless form is 

d_ 6p c (t) + 6p c (t) ■ «W Q (t - t q ) + 6W f (t - x f ) (57) 

dt 

+ n [ 5p (t) - 6p (t-T)] 
c c 

assuming the the pressure and temperature at any given instant are constant 
throughout the combustion chamber and the time lag is constant for all propellant 
elements. The dimensionless chamber press’ire, P Q and flow rates w Q and Wf are 
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defined in terms of their steady-state values and the dimensionless time, 
t, time laps for the oxidizer and the fuel, t q and t^, and sensitive time 

lag, ~, are defined in terms of the gas residence time, e • (Symbols are 

& 

listed at the front of the report and 6g is defined in Appendix A. ) 

Two dimensionless equations representing the dynamics of the feed 
system, one for the oxidizer and one for the fuel assuming constant feed 
pressure at the injector inlet and all capacitance and inertance located at 
the injector are 

2 2 

P o (l+J o ED ) 6p c (t) + (l+JD + J q E q D ) 6* o (t) = 0 (58) 

2 2 

P f (l+J f E f D ) 6p c (t) + (1+JjP + ) «W f (t) = 0 (59) 

where D is a differential operator, and the dimensionless inertance, J Q and 
Jf, capacitance E 0 and E f and pressure drop parameters, P 0 and Pf, are defined 
in Appendix A. Substituting equations 2 and 3 into equation 1 and applying 
Laplace transformation to the result yields 

2 2 
[S+l-n+ne' XS ][l+J S+J E S ] [ 1+JS+JES ] 

OOO I X I 

= -e 0 P o [l+J o E o S 2 ] [l+J f S+J f E f S 2 ] - e TfS p f [l+J f E f S 2 ] [l+^S+J^S 2 ] (60) 

Substituting a+iw for S in equation (6o) and equating the real and imaginary 
parts of both sides results in two simultaneous equations , 


-ot . 

(GH-KN ) ( o+l )- ( KH+GN )a> - -P e ° (coo ut [MH-N(2awJ E )] + sin arc [MH+H(2cuuJ E )]\ 

o I o oo o oo J 

« t % 


f f ) 

-P f e | coo uT f [RG-K(2ouJ f E f )] + sin wt f [ KR+G( 2a<uJ f E f ) ] j 


(61) 


NJh L M 

(GH-KN)w + y KH+GN ) (a+l) = -P Q e 0 J coo ut o [MH+H(2cuuJ E )] 

-sin ait [MH-N(2awJ E )] } 
o o o J 

-F f e j cos u)T f [KR+G(2oaxJ f E f )] -sin uix f [RG-K(2ouJ f E f )] J 


(62) 
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where : G 

H 
K 
N 
M 
R 


1 + J a + J E (a 2 - w 2 ) 
o o o 

1 + J f a + J f E f (a 2 - to 2 ) 

J to + 2 J E otto 
o oo 

J f a) + 2 J f E f ato 

1 +J E (a 2 -to 2 ) 

0 0 

l+J f E f (a 2 -u) 2 ) 


For reasons to be stated later, the interaction index n and the 
vaporization time t have been set to zero in Equations (6l) and (62). 

Equations (6l) and (62) can be solved simultaneouslyf for <* and w to 
evaluate the stability of any combustor design once the pressure drop, in- 
ertance, capacitance and time lags axe known*. The magnitude of a indicates 
incidentally, the proneness of the system to instability (if o>o) or to stability 
(if a<o). 

The time lags, t q and x f are by Crocco's definition the total time lags 
and are composed of a constant, steady-state (insensitive) portion and a var- 
iable (sensitive) portion x. Rigorous analyses would take into account this 
time -dependent or sensitive portion of the time lag; however, for most appli- 
cations it can be neglected because it is small compared to the total time lag. 
Once the sensitive time lag is taken as zero it follows that the interaction 
index must also be zero because zero sensitive time lag requires that the 
burning rate be independent of chamber pressure. 


t The programs written for this purpose are listed in Appendix C. The 
calculations which underlie the solutions for a and to are discussed in 
Appendix B. 

* The curves showing variations of a with L 0 (V 0 ) while Lp (V F ) is held 
constant are shown in Figures C (d). 
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The total time lag is defined as the time increment between injection of 
a propellant and its conversion into combustion products. Certainly this time 
lag is not the same for all propellan . elements. It is therefore customary to 
define an average time lag for each propellant which may be done by determining 
the lapsed travel time between injection and the axial position where 
combustion is assumed to take place. The methods of Priem [7] can be used 
to determine the position of the combustion front. 
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DISCUSSIONS AND RESULTS 

In to illustrate the effects of capacitance, inertance and 

resistance of ti.e systems on stability, equations 5 and C were solved for 
the neutral stability boundaries (ot = 0). Figures 1 through 3 show the 
neutral stability boundaries in terms of the fuel and oxidizer orifice lengths 
for three different oridizer pressure drops. Fuel and oxidizer cavity volume, 
fuel pressure drop and time lags are held constant. Figures 4 and 5 show 
the neutral stability boundaries for two different fuel pressure it' ops 
while oxidizer pressure drop is held constant along with the cavity volumes 
and time lags. Because more than one pair of roots satisfy equations 5 and 
6 multiple stable and unstable zones exist. The fuel and oxidizer pressure 
drops affect these zones but unfortunately nc trends are apparent. 

Figures 6 through 8 show the neutral stability boundaries in terms of 
the fuel and oxidizer cavity volumes for three different exidizer pressure 
drops. Fuel and oxidizer orifice lengths, fuel pressure drop and time lags 
are held constant. Figures ? and 10 show the neutral stability boundaries 
for two different fuel pressure drops while oxidizer pressure drop is held 
constant along with the orifice lengths and time lags. Although multiple 
statle and unstable zones exist as in figures 1 through 5, it appears that 
intermediate values of oxidizer and fuel pressure drops (figures 7 and 9 , 
respectively ) result in the narrowest unstable zones. Thus, if an operating 
point were selected which is in an unstable zone of figure 6, the system 
could be stabilized by reducing either the oxidizer pressure drop to 124 psi, 
or the fuel pressure drop to 87 psi (figures 7 and 9). Further reduction in 
either oxidizer or fuel pressure drop results in the return to unstable 
operation. This result cannot be predicted by either single time lag 
models or double time lag models which do not include injector inertance 


and capacitance. 
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CONCLUSIONS 


The contribution advanced by the stability model described here is 
the inclusion of injector inertance and capacitance in addition to resis- 
tance in the analysis of bipropellant rocket systems with different time 
lags. Neutral stability boundaries are shown in terms of these parameters 
in order to demonstrate their interactions. 

This model provides a method of designing a stable system by optimizing 
the pertinent design variables rather than maximizing the pressure drop and 
ignoring the others which has been the traditional approach. This model 
suggests that in some cases stability can be enhanced by reducing pressure 
drop, and therefore maintaining pressure drop to chamber pressure ratio is 
not necessarily desirable. 
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Subscript 


Combustion chamber 
Oxidizer 


Fuel 


Initial 


Superscript 


Indicates dimensional quantities 



LIST OF SYMBOIS 


SYMBOL 


LKFINITION 


P 

E 

J 

P 

P* 

P 

M 

V 

T 

• 

m 

R 

P 

P 

X 

e 

g 

1 

a 


s 


w 

u 

n 

a 

o» 

u 

w* 

T 

T* 

T 

T* 


Differential operator, d/ctt 
Dimensionless elasticity paLrameter 
Dimensionless inertia parameter 
Dimensionless pressure drop parameter 
Pressure 

Dimensionless pressure 

Molecular weight 

Volume 

Temperature 

Mass flow rate 

Universal gas constant 

Pressure drop 

Density 

Compliance 

Gas residence time 

Injector orifice length 

Injector orifice area 

Laplace transform variable 

Weight flow rate 

Dimensionless time 

Interaction index 

Dimensionless Damping Coefficient 

Damping Coefficient 

Dimensionless frequence 

Frequency 

Dimensionless time lag 
Time lag 

Dimensionless sensitive time lag 
Sensitive time lag 


Appendix A 

Definition of Parameters 


The gas residence time, C_ is used to nondimensionalize many of the 

& 

quantities used in the analysis. It is defined as 

V p V p M 
c c _ c c c 

g m = R T m 
c 


The pressure drop parameter, P is defined as 


The elasticity parameter, E is defined as 
E - 

E me 

g 


The inertia parameter, J is defined as 


2ApA6 


The read, and imaginary parts of the Laplace transform variable S (o and u>), 
the oxidizer and fuel combustion delays, and the sensitive time lag are non- 
dimensionalized by the use of the gas residence time, $g, namely, 
a = a*0 

g 

(jj ss (tf*0 


t = t /e 

0 0 g 


'f = V 9 e 

t = T*/e 


i * 



Appendix B 


Solution of Equations (6l ) and (62) by Nevton-Raphson Method 

For convenience , Equations (6l) and (62) are denoted, symbolically and 
respectively, as 

f(a,o>) = 0 (A-l) 

g(a,u>) - 0 (A-2) 

By series expansion. Equations (l) and (2) are rewritten 


3f 


= flvi-Vi 1 * (a-0 n-i > to ( VrVi 

+ (u-u ) 7— (o ,<i) ] 

n-1 3w n-1* n-1 


) 


(A-3) 


g(a,w) = g(a ,w ) + (a-a ) (o ,w ) 

n-1 n-1 n-1 3a> n-1 n-1 

+ (w-u , ) (a ,w ) 

n-1 3u n-1 n-1 


(A-l*) 


where « n _^, n _^ are the values of a and u at stage of computation n-1. 

Now if at a subsequent stage n, where a = a and w = u> the right-hand sides 

n n 

of (3) and (4) vanish, the following can he written immediately. 

°n to < Vl’Vl ) * “n to 1 “n-1 '"n-1 1 * - f( Vl>“n-l ) 

+a n-l to •Vl'Vl' 

%-i E <Vi'Vi> (A - 5) 


a T® - (o , ,w , ) + u (a . ,u „ ) » -g(a .a) 
n 3a n-1 n-1 n 3u n-1* n-1 ; ® ' n-1 n- 


1) 


*Vi to ( Vi'Vi> 


*Vi to ( Vi'Vi> 


(A-6) 



Written in matrix form, (S) and 


df 

af" 


° n 

ao 

3u) 


FI 

11 

31 


“n 

do 

d(i> 

n-1 



Inversion of (7) yields 


- - 
°n 


• 

af 

* 

af 



do 

du 


a _ 



un 


11 

11 

ft 4 


do 

dw 

1 


• 


take the form 



21 

do 

21 

dw 


Vl 

♦ 





n-1 

11 

11 


“n-1 

do 

dw 

• 

n-1 

» « 


■ « 

f 

♦ 

Vi" 

9 

> « 

n-1 

“h-1 

» 4 



1 

m 


21 

du 


where 


DFS • 


Vi 


L 

j 



21 

21 

c.« 

a«» 

21 

21 

do 

du 


_21 

da 


- 21 


* * 

f 

dtt 



?f 


9 

do 

m 

n-1 

. . 


(A-7) 


(A-8) 


(A-9) 


Iteration follows until the roots a and » are pin-pointed within desired 
accuracy. 

Such a scheme of solving simultaneous equations Is known as the Newton- 



Kaphson method. 

A computer program was developed such that for given parameters for 
a specific rocket engine, for instance, 


P 0 “ 

•‘V 

«< 4.82, 

J 0 - 1.189, 

% - 0.2018, 

t 0 - 5.269 2. 

T r 

* 0.4. 

J f - 0.061, 

B f - 0.233, 

the fundamental roots 

are 

obtained 




u » - 0.112 
u. 0.49. 

Other examples are solved similarly and the results of computation 


are shown in Figs, i and ii 










APPENDIX C 


Computer Programs and Sample Calculations 


T 





1 

? 

3 

4 
4 

ti 

1 

9 

Q 

10 

11 

12 

H 

14 

14 

14 

17 
I 1 
1 * 
27 

21 

22 


21 

24 

25 
2 * 
27 
2 ** 
29 
10 


♦ I 0 H 

f B|PR )P6ll ANT 

C NAkCM 29, 1973 

r 

r CURVE 3, FIG. 1, PER MARCH 9, DATA FROM RICHMOND 
C SEEKING A PLOT OF LO VERSUS IF 

N« 1 

''RArMStlOOl NC ASF 
109 FORMAT ( f 2 ) 

1 RfAOCSaOU RLO,«LF 
101 c flRMAT ( 2F 10* 6) 

r 

GC*32.2 
THG*0, 00026 

r 

OPO-62. 

AT0#0** 

RMO-25.* 

C 

DPF**3.3 

ATFnO.945 

*MF*4.30 

c 

poo*io.o8 

PFFsA.RO 

PJO»RLOPRMO/(24.*DPn*ArOACCATHGl 
PJFsRtFPRMF / ( 2 4 • *nPFAATp*GC*TMGI 
RFo#o*07? 1 
PFF* 3.241 
r 

TO*5.3AS 

TF=O.S39S 

f 

C THE FOLLOWING CAPOS CONTAIN INFO FOR LATER USE IN VARYING WO VS VVF 

C 

C PCC* 1 248 • 

C Wn*12.4 

r RHnp=o#04 l l 

f VVF*9 . c 

C HHFR-0.000239S 

r i*o*o.ot 

C RKF-0.680 

r pfo*MHnn*RKn*vvn/|Rwri9THGP|Pon*o*5| » 

r PFF*«HrF*RKF*VVF/|RMF4THG*(PFF40*5ll 

t POO.S.04 WAS USFf) IN CURVE 2 f FIG.l 

r OPO-124 WAS USFD I IN CURVE 2« FIG.l 

f *)PF *94,6 WAS USED UN CURVE 2, FIG.l 

r. Pf 0*0* 1 3R 7 WAS USEO IN CURVE 2'FIG*1 
r FNO OF LATENT IF OR MAT I ON 

r 

Al =0*01 
W*O.Ol 

r 

17 FORMAT C IH 1 1 

WRt TE I 4* 1 7) 

19 FORMAT C9X * •LINEAR ENGINE COMBUSTION STABILITY STUDY*,/! 

WRITF(6ilR) 

20 FORMAT(iOX, *Ln«*,El0.4 f SX,*LF«* ,E10.4,/) 

WRITE! 4,20) RLO.PLF 




M 1 PM R H 4 rU?X,ML*»HA' f llX»'ALPMA l • • 15X , • 3HFGA 1 1 10X» • OHEGA t • ,✓ > 

*? tTMTEfAfl?) 

33 S FT^l. 

34 «-V"‘A = rxP(-r 3*AU 

3* fvfa-^fxpi -rr*au 

34 i»pi: r »erosc n*wi 

)> PPSa*SIN(TO*w) 

3 ’ ‘’oCf =CnS( TF*W) 

3Q PPSF*S!NI TF*W) 

40 PPC T-* I • 

41 PPST«0. 

4? A AWW-Al *AL-M*W 

C 

43 PGM .♦Pjn*4l>Pjn*PE0*AAWW 

44 PMsl.4PJF*AL*PJF4PEF*AAWM 

45 PK*PJ04W4 Cl. ♦?.4P€0*AL» 

46 *NN*PJF*WM I • *2 • *PEFPAL I 

47 PM»or,-PJO*AL 

♦ 4 PP*P4-PJFPAL 

4-3 PU=2.*Pjn*PEn*AL*W 

50 PV*2 • 4PJF *PFF*AL*W 

41 PA»AL*U-PN*PN*FTA4PPCT 

52 PR*W-PN*F TAPPPST 

43 CC = FVOA*PPCn 

54 FFfcFVOAPPPSO 

55 1>0=CVFA*PPCF 

56 FF-*EVF A4PPJF 

57 AJl*»G*PH-PNN*PK 

55 AJ2 3 PH*PK^PNN*PG 

50 A J5«PH*PH-PU*PN f < 

60 AJ7«PM*PNN*PH*PU 

41 AJQ=PR4PG-P^*PV 

6? A J 1 1 *PP*PK4PG*PV 

43 PGA*Pjn*< l.f2.*PFH6AU 

64 PMA*“PJF*( 1 • f2« *PFFPAL I 

44 "KA*2.*PJ0*P c n*W 

66 PNA«2* 4p JFPPEF *W 

47 °MA*2 • •PJ0*PF r l4At 

4« PQA*?.4PJF*PFFMl 

M PUA=PKA 

70 PVi'PMA 

71 PA A* 1 • 

72 PBA* 0* 

n cca*- Tn*fvoA*pocn 

74 FE A«- TO*FVOA*PPSO 

75 PDA«- TF4EVFA4PPCF 

76 FFA«- TF*€VFA*PP$F 

H PGH^-PKA 

’F PHW*-PNA 

79 PKW-PGA 

»0 PNW-PHA 

41 PMy**r>GW 

8? PRW-PHK 

83 PUW»P#A 

«4 PVK»PRA 

55 PAW«-P)3A 

86 PftW*PAA 

57 CCH-FEA 

59 EEW—CCA 

89 0 nw»FFA 




~>0 

r 

‘74 

n«j 

97 
on 
Q9 
1 )3 
1 91 
! 0> 

1 Oi 

104 

I OS 


c 

IQS 

C 

r 


F Fw--I)*)A 

CK l -f>G*PHA*PGA*PH-PNN*PKA-PNA*PK 
rK^=PH*PKW+PHW*PK+PNN*PGW>PNWPPG 
CK J*PG*PHW4-PGW*PM-P*N*PKW-PNH*PK 
rK4=PH4PKA»PHA*PK»P^N^PGA+PNA*PG 
CKS="HA*PHiP**®HA PIJ*PNA-PUA*PNN 
CK4* JM W *ph+ph*Pss-PUW*PNN-PU*PNW 
CK7sPMA+ 0 NN+PM*P^A>PH*PUA*PHAPPU 
CK A=PMrf*PMN+PM*PNW+PHW*PU+PH*PUW 
C*9=®R4*®G+PR*PGA-PKA*PV-PK*PVA 
CK10 = PPW*PG«-PR*PGW-PKW*PV-PKPPyW 
CKU *PRA*PK«-P«*PKA*PGA*PV+PG*PVA 
CK!2=PRW*PKfPR*PKWfOGW*PV*PG*PVW 

FFF = PA*AJl-Pa*AJ?*Pnn*tCC*AJ5*Ee*4J7UPFF*<0D*AJ9«-FF*AJill 
GGG=PR*AJ l*PA*A |?*P00*( CC*A J7-EE*AJ5 l-**PFFP( OD*AJ l l-FF*AJ9l 
F' FA=*PAA* AJl + PA+CK t'-PBA*AJ2-PB*CK4 
1 fpnn* ( CCA*AJ5*CC*CKS4-F6A*A J7*EE*CK7) 

1 fPFF*(00A*AJO400*CK9*FFA*AJU+FF*CKlt| 

FNO OF FFF A 

FFFW*PAW*4JUPA*CK3-P0W*AJ2-Pfl*CK2 
l ♦PnO*lCrw*AJ5*CC4CK6*EFW*AJ7*EE*CK8) 
l *PFFP< r)QW*AJ9^0U*CKl0>FFW*AJU4FF*CKl2) 

FNn nF FFFH 


1 0 7 GGr,A«P8A*AJUP9*CKUPAA*AJ2*PA*CK4 

L ♦POO*(CCAPAJ7»CC*CK7-ffa*AJ 5-EE*CK5I 
1 ♦PFF*( DDA*A Jt 1 ♦00*CK1 1-FPA*AJ9— FF*CK9 1 

C FNI) HF GOO A 
C 

l OP GGGW=PRW* AJ l ♦P8*CK3*-P AW*A J2*PA*CK2 

1 ♦P0n*ICCM#AJ74CC*CK8-EEW*AJ5-EE*CK6l 
1 ♦PFF*|OOM*AJHFOO*CK12-FFW*AJ9-FF*CKtO) 
r CNO flF GGGH 

r. 

100 OFG«FFFA*G^Gw *GGGA*FFFW 

110 GGG1*(GG ’FF-FFFW^GGG) /dfg 

111 A l* AL-GF' 

112 GGG2M-1 .GAAfFF4-FFFA*GGG)/0FG 

113 W1-W-GGG2 

r 

114 lb FilRHAT(10X,2E16.rt*5X, 2E16.8) 

US W»lTF(Atl6> AL * 41 » W*W1 

114 f F(ARS(GGGl > «LT. ( 0*001 )«AN0«ABSfGGG2)*LT* (0*001 1 1 GO TO SO 

IU Al*Al 

UP w » wi 

l lo G° TO A 

12) so If (NCASF-N)52»52 f 51 

121 SI M*Nf l 

122 Gn TO i 

12) S2 WR1TF (* f m 

124 Stop 

1 ** r N0 


if N try 



I ! if A« TNiWNf CUMUliS* l-IN STABILITY STUDY 


LUO. 10006 01 
ALPHA 

0.10000000F-01 
0.5L556450E 00 
-0.4*484670F 00 
-0.30994270F 00 
-0. 10960* 70F 00 
0.15632H70E 00 
0.89H94940E-01 
0.4H7780AOE 00 
-0.*'7U2‘>49«: 00 
-0.P232I40F 00 
-0.7A99SHQOF-02 
-0. 1 lb 1 1030E-03 
-0.477H940F-02 


lf-o.soooe 00 

ALPHA 1 

0.51556450F 00 
-0.494846706 00 
-0,i099*270e 00 
-0* 10960670E 00 
0.15632870F 00 
0.89A989406-01 
0.497700406 00 
-0.27U25406 00 
-0.U232140F 00 
-0. 78895090E-02 
-0.11*110306-03 
-0*477549406-02 
-0.473056706-02 


OMEGA 

0.100000006-01 
0. 417432 70E-01 
0.6983953QE-01 
0*729697306-01 
0*895799906-01 
0*157071806 00 
0*865384606 00 
0*376352606 00 
0*310356106 00 
0.33325400E 00 
0.390476106 00 
0.436219506 00 
0*436357406 00 


OMEGA 1 

0*417432706-01 
0.698395306-01 
0.729697306-01 
0.855799906-01 
0.157071806 00 
0.865384606 00 
0.376352606 00 
0.31035610E 00 
0.33J25400F 00 
0.39047810E 00 
0.4362 1950E 00 
0.43635740E 00 
0.43636080E 00 



pLta 




* jup 

r M | pk l n FL LAN F 

C THIS IS CURVE 2, FIS. 4, PER DATA FROM RICNMOND ON MARCH 8,1973 

t„ SrFKlNf, A PLOT OF VO VFRSUS VF 

I N = L 

» PF AO (5*100) NCASC 


3 

100 

FORMAT (12) 

4 

1 

RFA0C5, toll VV°,VVF 

r, 

101 

FORMAT (?F10*4) 

JS 


PN-O* 

f 


GC»12.? 

<* 


THG»0. 00026 

P 


ATO*0.6 

10 


ATF«0* 945 

11 


9X0*0.01 

12 


PKF»0*680 

11 


RMO*?5*4 

1'* 


RMF=4.30 ■ 

15 


c*,H0O*0.04ll 

16 


RHFT = 0*0002185 

17 


npn*i24. 

IR 


0PF*H6.6 

to 


pqd*i.16 

20 


PFPa 7.21 

21 


P JO* 1.1 09 

22 


PJF=0.747 

21 


PFO=RHnn*RKO*VVO/t RMO*THG*( P00*0* 51 ) 

24 


PFF=RHFF*RKF*VVF/(RMF*THG*(PFF40.5M 

25 


T0*5. 385 

26 


TFj*0 . 5 305 

27 


AL*0*01 

?R 


WaO.Ol 

29 

17 

FORMAT ( IHl ) 

10 


WR!TF(6,17> 

It 

18 

FORMAT 1 9X , 'L INEAR ENGINE COMBUSTION STABILITY STUDY', /) 

12 


WRTTF(6,1B) 
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